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An impurity solver for the dynamical mean field (DMFT) study of the Mott insulators is proposed, 
which is based on the second order perturbation of the hybridization function. After carefully 
benchmarking it with Quantum Monte Carlo results on the anti-ferromagnetic phase of the Hubbard 
model, we conclude that this impurity solver can capture the main physical features in the strong 
coupling regime and can be a very useful tool for the LDA+DMFT studies of the Mott insulators 
■ with long range order. 
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> ! I. INTRODUCTION 

o . 

Mott insulators are insulators which are caused by electron-electron interaction rather than the crystal structure, 
t — ' In sufficient low temperature, long range spin or orbital order will be established to get rid of the remaining entropy 
induced by the unfrozen spin or orbital degree of freedoms. Various types of interesting long range order may appear 
t—( ' in different systems, i.e. the anti-ferromagnetic (AF) order, the non-collinear spin order as well as the orbital ordering 
^ , states. The fruitfulness of the long range orders in the Mott insulator phases is a consequence of the detailed balance 
Ch ■ among electron kinetic energy, interaction energy and the spin orbital coupling in the solid. Therefore the ab initio 
c/2 \ study of the Mott insulators is very crucial to achieve quantitative understanding about these materials. 
,i_S ■ While the first principle calculations for the Mott insulators based on local density approximation (LDA) always 
encounter problems generated by the strong correlation effect among those partially occupied localized orbitals in 
these materials. For example, the LDA calculation incorrectly predict two typical Mott Insulators NiO and FeO to 
be metallic. This is due to the non-interacting form of the LDA type wave function can not capture the electron 
correlation effect, especially the strong on-site correlation effect in the relatively local orbitals such as 3d or 4/ orbitals. 
Although LDA underestimates the strong correlation effects in these materials, it still describes quite well the long 
\ range Coulomb interaction and those less correlated s p bands. Therefore a very efficient way to solve this problem 
- 1 ' is to implement the LDA with some many body techniques which can treat the correlation effect better. The first 
example of such scheme is the LDA+U method jl;], in which the Hartree-Fock mean field treatment is adopted to 
capture the orbital physics. The LDA+U method can correctly describe many Mott insulators, i.e. the YVO3 and 
' LaVOz 0- While since the wave function in LDA+U is still the single slater determinant, it can only capture the 
static orbital or spin correlation but not the dynamical correlation effect. Another attempt is the LDA+Gutzwillcr 
method proposed very recently by two of the authors, where the Gutzwiller variational approach [3[ Q is used to take 
into account the correlation effect (LDA+G). LDA+ Gutzwiller is superior to LDA+U in the sense that the dynamical 
7-H ■ correlation effect can be taken into account by the multi-configuration nature of the Gutzwiller trial wave function. 
t-H ' However, LDA+G is still very difficult to be used in the finite temperate problems because it is a variational method. 

In the past twenty years, the dynamical mean field theory (DMFT) Q has been quickly developed to be a powerful 
method to solve the strongly correlated models on the lattice. DMFT maps the lattice models to the corresponding 
^ . quantum impurity models with the environment to be determined self-consistently by the DFMT self consistent 
equations. DMFT keeps the full local dynamics induced by the local interaction and has been successfully applied to 
various of correlation systems, such as the Mott transition in Hubbard model @ @, the pseudo gap behavior in high 
Tc cuperates H as well as the heavy fermion system Q Combined with LDA, the LDA+DMFT has quickly 

become a very useful numerical tool for the first principle studies of strongly correlated materials, which has been 
successfully applied to many strongly correlated materials. 

The nuclei of DMFT is to solve the Anderson impurity model introduced by the DMFT self consistent procedure. 
The solver could be either numerical and analytical. The most popular two numerical solvers are exact diagonalization 
[ill ] [l2| method and Hirsh-Fye quantum Monte Carlo [l3| [3] [HI- Besides, the continuous-time quantum Monte 
Carlo (CTQMC) solver [l7| [H[ [HI] [13] has also been developed recently and greatly boosted the progress in this 
field . However, the numerical methods are often very time-consuming and for large systems it is better to use 
some analytical solver which is much faster and can capture the basic physics as well. The analytical solvers include 
Hubbard-I approximation [2l[ [22j , the equation of motion (EOM) method [23[ , the iterative perturbation theory [24j 
[IH, the fluctuation exchange approximation(FLEX) [2(| and the non-crossing approximation (2?| [28J. Among them 
the Hubbard-I solver, which simply takes the atomic self energy, is very convenient and easy to be implemented in 
the existing ab initio codes. Since it can capture the atomic multiplet effect (Hubbard bands) in a very efficient way, 
the Hubbard-I solver has been widely used to study the paramagnetic Mott insulators. While the Hubbard-I solver 
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can not be directly applied to the Mott insulator states with long range order, simply because the atomic self energy 
does not contain any information about the long range order in the system. The way to solve this problem is to go 
beyond the atomic approximation and include the leading order corrections from the heat bath as well. 

A second order strong coupling expansion solver has been developed by one of the authors in reference [29[ . Although 
the long range ordering temperature can be correctly calculated using this solver, it also suffers from the following 
facts, i) Because the strong coupling expansion is based on the atomic limit without long range order, it is not valid 
to treat the systems with large order parameter, i.e. the AF phase in low temperature with the fully developed AF 
order. ii)Due to the multiple summation over the configuration space, the calculation is still quite heavy for the 
multi-orbial system especially the /-electron system. 

In order to overcome the above two problems, in this paper we propose a new scheme to include leading order 
corrections from the heat bath. The idea is quite simple: we use the second order perturbation theory to correct the 
energy level of each atomic eigenstates. Then we construct the "atomic" Green's function and the corresponding self 
energy on a single site based on the re-normalized atomic eigen-energy. With the above procedure, we can obtain 
the main contribution of the heat bath to the self energy in the Mott insulator case, which is the re-normalization 
of the eigen-energy of the atomic configurations. Because only the eigen-energy not the occupation of the atomic 
configurations is treated by second order perturbation, the present method can be applied to very low temperature 
with quite large order parameters. With the present impurity solver the static inter-site spin or orbital correlation, 
such as the super-exchange process, can be taken into account. Compared with the previous atomic expansion solver, 
this solver does not have multiple summations over the atomic configurations, which make it very fast and is more 
suitable for the study of the /-electron systems. 

The paper is organized as the following. We will present the detailed formalism of the impurity solver in section II. 
After that in section III, we will benchmark it on both paramagnetic and anti-ferromagnetic Mott insulator phases 
in the Hubbard model. And finally we will carefully discuss the results and make the conclusion. 



II. DERIVATIONS 



A. General Formalism 



The crucial part of the dynamical mean field theory is to solve the quantum impurity model with the following 
Hamiltonian: 
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The above expression is a general form of multi-orbital Anderson impurity model with arbitrary local interactions as 
well as the heat bath introduced by DMFT self consistent procedure. In this Hamiltonian, m, m! and n, p, q are the 
combined spin and orbital indices, while / denotes the energy levels in the heat bath. 

In the LDA+DMFT study, we often need a fast impurity solver which can capture the most important features of 
the impurity problem. For the Mott insulators, such a fast solver can be constructed by doing perturbation around 
the atomic limit. In the present paper, we apply the "contour integral" technique to obtain the electron Green's 
function [27j ■ The retarded form of the electron Green's function on the impurity site can be written as 

GY n \u +i„) = i § c fee-** UF m U(F"^)p a 

a/3 

xPp(z)Hp^ a (z, Z + UJ + \T])P a {z + UJ + 1VI) (2) 

Z = § D &e- f> 'Y.P*{*) (3) 

a 

In the above equations, c is a contour in the complex plane, oriented counterclockwise, and surrounding all singularities 
of the resolvent (z — Hi oc )~ l , where Hi oc is the local part of Hi„ lp ; a denotes the ath eigen state of Hi oc with the 
eigen value of E a ; (F m ) a p is the matrix element of local fermion operator f m ; Hp a (z, z + u> + 177) is a proper vertex 
function to be discussed below; Z is the partition function of the impurity problem and P a {z) can be viewed as the 
"propagator" of the atomic eigen states. In the absence of hybridization takes the simple form of P a ,o( z ) — ^/( z ~E a ). 
When taking the hybridization part into account, P a (z) can be expressed as 

P a (z) = l/(z-E a -X a (z)) (4) 
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FIG. 1: The diagrams of the lowest order contribution to the self energy of local atomic states. The break lines denote the 
propagators of atomic states, the continuous lines the propagator of electrons in the heat bath, and each vertex corresponds to 
an element of a local Fermion operator with the coupling strength V . 



where S Q (z) can be viewed as the self-energy of the cvth local configuration that comes from the hybridization between 
the impurity site and the bath. 

In order to evaluate E tt (z), we use the method heuristic from the techniques in the large-N expansion. The self 
energy can be expressed in a diagrammatic way as it is introduced in detail in the review [30(. In this paper, we focus 
on the Mott insulator system, which only requires the lowest order diagram as shown in Fig [I] 

The expression of self-energy corresponding to the above diagram is: 

p i mm' ^ 
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p I mm 
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where /(ez) is the fermi distribution function, while 



Aw(^ + i7?) = V V <™ Vlm ' (6) 
lo + ir] - ei 

is called the hybridization function in the dynamical mean held theory and its imaginary part is denoted as 

A^ m ,(e)=Im[A w (e + i77)] (7) 

in equation ([5]). 

Afterwards, for the sake of a simplest self-consistent approximation that provides a soluble set of equations, the 
vertex function in Eqn.@ could be set 

n a , a ,{z,z') = l (8) 
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and we now get the final expression of G™ m (u + irj) in Eqn.([2|). 

The local electrons' quasi-particle spectrum functions are the direct result of the retarded form of the Green's 
functions. For a given Hamiltonian as Eqn.([T]), one may use the equations above to evaluate the Green's function. 
Implemented with the DMFT self-consistent condition described in the review this method can be used as an 
impurity solver in the DMFT study of the multi-orbital systems with complicated local interactions. 



B. Single Pole Approximation 

The method shown above works well for Mott insulator phase, and a few results will be shown in the next section. 
However, here we introduce a further approximation which makes the impurity solver much faster and more convenient. 
The main idea of this approximation can be interpreted as the following. For a Mott insulator solution in DMFT 
with long range order , the most important effect of the heat bath is to modify the effective energy level of the atomic 
configurations. Therefore in the present paper, we only consider the effect of heat bath in the renormalization of the 
atomic levels, but neglect the possible broadening of the atomic levels. 

With the above consideration, Eqn.Q can be approximated as 

P a (z) = l/(z-E a -E a (z)) 

« l/(z-E a -Re[X a (E a )]) 

= l/(z-E a ) (9) 

, and E a = E a + He[T, a (E a )] is the renormalized atomic level containing the second order correction of the hybridiza- 
tion function. Because the Green's function P a (z) contains only one pole at E a , as the same form in unperturbed 
case, we call it "single pole approximation" . Thus the electron Green's function can be written in a similar form as 
the atomic Green's function 

G? m '(w + ir,) =^E(i ?m W(Wt)^ (10) 

' a/3 

where the partition function is Zf = J2 e ~ ■ Correspondingly, the self energy that comes from local two-particle 

a 

interactions becomes 

E(w + irj) = oj + ir) + (i - 1 - Gj 1 (12) 

This self energy can be implemented in the DMFT self consistent loop and we obtained a fast impurity solver for the 
Mott insulators with the long range order fully developed. With the current approximation, the DMFT here is equiv- 
alent to a generalized static mean field approximation using the local many-body bases. The spontaneous symmetry 
breaking can be well described by the difference in energy and occupation for atomic configurations which can be 
connected by symmetry, i.e. the singly occupied spin up and down states in the one-band Hubbard model. Therefore 
the super-exchange process can be thus captured through the self consistently determination of the hybridization 
function in the DMFT loop. 



III. RESULTS AND DISCUSSION 

A. Single band, Paramagnetic phase 

The spectral function of one-band Hubbard model with various U at half filling is plotted in Fig [2j It works on the 
Bethe lattice with the semicircular density of states 



D{e) = -^y/U*-#, \e\ < 2t (13) 

where t = 0.5 as usual. Here we do not use the "single pole approximation", which means Eqn.Q is implemented. 
As a perturbation method around atomic limit, we echo the fact that this kind of solver is limited to integer filling 
and large U [2^|, that is, the present solver is just good to study in the regime of Mott insulating state. 

As a comparison, we also exhibit the result of present solver with and without single pole approximation compared 
with that of DMFT+NRG, see Fig [3] It shows that the result of present solver at large U regime agrees well with 
that of NRG, and the single pole approximation captures the main features of the Hubbard bands. 



FIG. 2: The spectral function of single band Hubbard model on the Bethe lattice at half filling for various values of U, the 
energy unit is chosen as D=2t=l. 




FIG. 3: The spectral function for U = 4 on the Bethe lattice compared with DMFT+NRG, the energy unit is chosen as 
D=2t=l. 



B. Single band, Antiferromagnetic-Paramagnetic Phase Transition 

In this section, we study the antifcrromagnetic(AF)-paramagnctic phase transition in single band Hubbard model 
at half filling. We focus on the AF order parameter m = (n^ — nj,) on a given sublattice the Neel temperature. 
Henceforth we use the single pole approximation, with the self energy in Eqn. (|12p . First we study the single band 
Hubbard model on the Bethe lattice with the semicircular density of state. Fig[4] shows the magnetisation m on one 
sublattice versus temperature T for U = 4 and U = 7. The curves show that the phase transition is of second order, 
and the critical temperatures are about T c = 0.061 for U — 4 and T c = 0.037 for U = 7 respectively. Likewise, we get 
T c under different U and plot them in Figj5] This result is comparable to that of DMFT+QMC in large U regime, as 
it is expected. 

We also apply this method to estimate the Neel temperature on the 3D isotropic simple cubic lattice. As the 
dispersion relation is E(k) = 2t(cosk x + cosk y + cosk z ) , we choose t = 1/6 to ensure that the half band width 
equals to 1. It is known that, in the large U limit, on the 3D cubic lattice, the Hubbard model is equivalent to a 
Heisenberg model with an antiferromagnetic coupling J, and the relation between the parameters reads J = 4i 2 /U 
[3l| . Accordingly, the correct Neel temperature of the model is Tjv = 3.83i 2 /i7 [32j while the Weiss mean field gives 
the result Tjv = 6t 2 /U 33]. Remember that we only take the lowest order of corrections to the local interactions, 
so the result should be close to the one of Weiss mean field, see Fig [6] We have also plotted in this diagram the 
Hartree-Fock Neel temperature for the same model (34|, which is completely wrong in the large U limit. 

It is intuitive to compare our results to that of Heisenberg model, because the latter has been studied intensively. 
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FIG. 4: The sublattice magnetisation m versus temperature T for U = 4, 7, the energy unit is chosen as D=2t=l. 
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FIG. 5: The comparison of the Neel temperature obtained by DMFT+QMC (thick line) and present solver (circle dots) on the 
Bethe lattice at half filling. The half band width is used as the energy unit by defining D = 2t — 1. 



For example, we have compared the antiferromagnetic condensed energy. In the antiferromagnetic spin-wave theory 
for Heisenberg model, this condensed energy is 12.96t 2 /U on 3D simple cubic lattice [35|, about twice as much as 
Tjv of Weiss mean field. Here we have calculated the energy difference between antiferromagnetic and paramagnetic 
solutions, at a temperature of 0.002r; which is sufficiently low and plotted it in Fig [7] together with the condensation 
energy obtained by spin wave theory for the corresponding Heisenberg model. Our results are in good agreement with 
the spin wave theory which indicate that our method works very well even in low temperature. 
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FIG. 6: The comparison of the Neel temperature obtained by Weiss MF and present solver on the 3D cubic lattice at half 
filling. For cubic lattice the unit energy is chosen as the half band width D = 6t — 1. 




FIG. 7: The antiferromagnetic condensed energy compared with the spin-wave theory for Heisenberg model, where the unit 
energy is D — 6t — 1. 



C. Two band, Antiferromagnetic-Paramagnetic Phase Transition 

We have also studied antiferromagnetic-paramagnetic phase transition in two-band Hubbard model at half filling, 
using DMFT with single pole approximation. Because there are two band on one single site, the situation becomes 
much more complicated. Here we only consider the simplest case. First, we consider two degenerate band, with 
semicircular density of states and the band width is set to be 1. Second, we neglect the hopping between different 
orbits. Therefore the Hamiltonian can be written as [36| 

H at — U 2J n b ^n b , i + U' 2J m,o-«2,o-' - J fix^h^a 

b <j,a' a 

(14) 

There are three adjustable parameters above, which are U the intra-band Coulomb repulsive energy, U' the inter-band 
Coulomb energy, and J the interaction between different bands with the same spin. For J > 0, the third term on the 
right hand side in Eqn. (fT4")) contains the Hund's coupling —J'SfS% since S z = l/2(n-f — njj. The discussions below 
are based on this Hamiltonian. 

Here we further assume U — U' = 2 J, which is true for systems with the cubic symmetry [371 ] . The Neel temperatures 
versus various U obtained by us for J = and J — 1 are shown in Fig[5] respectively. 

The Neel temperature of two band Hubbard model with U = U' is lower than that of one band model with the same 
U. This can be understood in the way that in a two band system at half filling, there is AF order developed in two 
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FIG. 8: Two band model on the Bethe lattice under the condition of U — U' = 2 J for J = and J = 1 respectively, compared 
with one band result. The energy unit here for the Bethe lattice is chosen as D = 2t — 1. 



different orbitals. But there is no term to lock the relative fluctuation between them, which reduces its tendency to 
be anti-ferromagnetic. With non-zero Hund's rule coupling J, such relative fluctuation between AF order in different 
orbitals will be greatly suppressed, which result in the dramatic increment of the Neel temperature as shown in Figj8] 



IV. CONCLUSIONS 



This paper presents a new impurity solver based on the perturbation theory in the hybridization around the atomic 
limit. For the sake of speed and convenience, we have introduced a further "single pole approximation" with which 
we carefully studied the antiferromagnetic-paramagnetic phase transition in one band and two band Hubbard model 
at half filling. The DMFT with the current impurity solver thus provide a way to do mean field approximation with 
the atomic configurations rather than energy levels for the single particle states. Comparing with the DMFT+QMC 
on the single band Hubbard model, the Neel temperature obtained by this new impurity solver fits quite well with 
the DMFT+QMC results, which indicate that is can be very useful in the LDA+DMFT studies for Mott insulators. 
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